package it.crosato.stage.server.model.phylogeny;

import java.util.Vector;

public class NJ {
	int N;			// Numero inziale di fogle
	int K;			// Numero di cluster
	ClusterNJ[] cluster;		// I cluster dell'albero risultante
	double[] r;			// La distanza media dalle altre foglie

	/**
	 * Costruisce l'albero a partire dalla lista degli organismi e dalla
	 * matrice di distanza
	 * @param ds matrice di distanza
	 * @param org lista degli organismi
	 */
	public NJ(double[][] ds, Vector<String> org) {
		N = ds.length;
		cluster = new ClusterNJ[2*N-1];
		for (int i=0; i<N; i++)
			cluster[i] = new ClusterNJ(i, ds[i], org.get(i));
		r = new double[2*N-1];
		K = N;
		while (K < 2*N-2)
			findAndJoin();
		int K2 = K-2;
		while (!cluster[K2].live())
			K2--;
		double dij = d(K2, K-1) / 2;
		cluster[K] = new ClusterNJ(K, cluster[K2], dij, cluster[K-1], dij, null);
		K++;
	}

	/**
	 * Calcola la matrice di correzione
	 */
	void computeR() {
		for (int i=0; i<K; i++)
			if (cluster[i].live()) {
				double sum = 0;
				for (int k=0; k<K; k++)
					if (cluster[k].live() && k != i)
						sum += d(i, k);
				int L = 2 * N - K;	// Il numero corrente di foglie
				r[i] = sum / (L - 2);	
			}
	}

	/**
	 * Ritorna la radice dell'albero
	 * @return il cluster radice
	 */
	public ClusterNJ getRoot(){
		return cluster[K-1];
	}

	/**
	 * Ritorna la distanza tra due cluster
	 * @param i indice primo cluster
	 * @param j indice secondo cluster
	 * @return distanza tra i cluster
	 */
	public double d(int i, int j) {
		return cluster[Math.max(i, j)].dmat[Math.min(i, j)];
	}

	/**
	 * Cerca la coppia di cluster da unire e li unisce
	 */
	void findAndJoin() { // Trova i due cluster live più vicini e li unisce
		computeR();
		int mini = -1, minj = -1;
		double mind = Double.POSITIVE_INFINITY;
		for (int i=0; i<K; i++)
			if (cluster[i].live())
				for (int j=0; j<i; j++)
					if (cluster[j].live()) {
						double d = d(i, j) - (r[i] + r[j]);
						if (d < mind) {
							mind = d;
							mini = i;
							minj = j;
						}
					}
		join(mini, minj);
	}

	/**
	 * Unisce la coppia di cluster in uno unico rimuovendo ed elimina i due cluster
	 * @param i indice primo cluster
	 * @param j indice secondo cluster
	 */
	public void join(int i, int j) { // Unisce i e j per formare il nodo K
		double[] dmat = new double[K];
		double dij = d(i, j);
		for (int m=0; m<K; m++)
			if (cluster[m].live() && m != i && m != j)
				dmat[m] = (d(i, m) + d(j, m) - dij) / 2;
		double dik = (dij + r[i] - r[j]) / 2;
		double djk = dij - dik;
		cluster[K] = new ClusterNJ(K, cluster[i], dik, cluster[j], djk, dmat);
		cluster[i].kill();
		cluster[j].kill();
		K++;
	}

	/**
	 * Controlla che l'albero rispetti la proprietà dell'additività
	 * @return true se rispetta tale proprietà, false altrimenti
	 */
	public boolean checkAdditivity() {
		for (int i=0; i<N; i++)
			for (int j=i+1; j<N; j++)
				for (int k=j+1; k<N; k++)
					for (int m=k+1; m<N; m++) {
						double dijdkm = d(i, j) + d(k, m);
						double dikdjm = d(i, k) + d(j, m);
						double dimdjk = d(i, m) + d(j, k);
						if (!(dijdkm == dikdjm && dijdkm >= dimdjk
								|| dijdkm == dimdjk && dijdkm >= dikdjm
								|| dikdjm == dimdjk && dikdjm >= dijdkm)) {
							return false;
						}
					}
		return true;
	}
}